#####
# This R code replicates the analysis in 
# Ahlquist, John S. ``Policy by Contract: electoral cycles, 
# parties and social pacts, 1974-2000.'' \emph{Journal of Politics}.
# Last updated 10/31/2009
#####

#required libraries
library(survival)
library(foreign)
library(MASS)


#save data to local drive and set local working directory

pact.data<-read.csv("AhlquistJoP2010RepData.csv")

pact.data<-subset(pact.data, year<2000)


pact.data$exec.party<-as.factor(pact.data$exec.party)
pact.data$exec.party<-relevel(pact.data$exec.party, ref="R")

pact.data.discont.risk<-subset(pact.data, under.pact==0) # identifies risk set
attach(pact.data.discont.risk)

## fact referenced in paper:
pact.under<-subset(pact.data, under.pact==1)
pact.not.under<-subset(pact.data, under.pact==0)
mean(pact.under$leftgs, na.rm=T)
mean(pact.not.under$leftgs, na.rm=T)


## Figure 1

cox1<-coxph(Surv(start,stop,status.pact.new)~1)
cox2<-coxph(Surv(start,stop,status.pact.new)~strata(maastricht))

surv1<-survfit(cox1)
surv2<-survfit(cox2)
years<-c("1975:I", "1980:I", "1985:I", "1990:I", "1995:I", "2000:I")
poss<-c(5,25,45,65,85,105)

par(mfrow=c(2,1))
plot(surv1,ylab="Survival", conf.int=F,xlab="year", axes=F,lwd=2,xlim=c(0,120),
    main="Expected survival and New Pacts")
axis(1,at=poss,labels=years)
axis(2)

points(surv1$time,rep(0,length(surv1$time)),pch="|")
text(nn$period,.05,srt=45,labels=nn$country, pos=4,cex=.85,font=1)

plot(surv2,ylab="Survival", xlab="year",
    main="Stratified by Maastricht", lty=c(1,3),xlim=c(0,120),lwd=2,axes=F)
axis(1,at=poss,labels=years)
axis(2)
dev.off()


## models for table 2
fit.1<- coxph(formula = Surv(start, stop, status.pact.new) ~ lag.inflation + lag.unemployment + lag.deficit + lag.curacct + grgdpch + lag.openk + strata(maastricht) + tsle.ciep + leftgs + log(enpp)+ as.factor(cab.maj)+ frailty(wbcode,dist="gaussian", method="aic",caic=T),x=T)

fit.2<-update(fit.1,.~. - leftgs + exec.R)

fit.3<-update(fit.1, .~. - tsle.ciep + elec.in.next.6q)


## figure 2

###Interpetation plot
## plot interpretation
# 1 additional qtr. is on avg 1/16 of a CIEP, or
# median qtrly change in unemployment is 0.20
# median qtrly change in inflation is 0.3
# median election-to-election swing in left seats 4.3

b.hat<-fit.1$coef
Sigma<-fit.1$var2
sim.b<-mvrnorm(1000,mu=b.hat,Sigma=Sigma)
bhat<-subset(sim.b,select=c(tsle.ciep,leftgs,lag.unemployment,lag.inflation))
bb<-matrix(NA,ncol=ncol(bhat),nrow=nrow(bhat))
xsim<-c(100/16,4.3,.2,.3)
for(i in 1:ncol(bhat)){
  bb[,i]<-bhat[,i]*xsim[i]
}



bb.bar<-apply(bb,2,median)
bb.upper<-apply(bb,2,quantile,.975)
bb.lower<-apply(bb,2,quantile,.025)
b.exp<-100*(exp(bb.bar)-1)
upper<-100*(exp(bb.upper)-1)
lower<-100*(exp(bb.lower)-1)

yrng<-c(100,75,50,25)

pdf("coxinterpJoP.pdf")
plot(b.exp,yrng,ylim=c(10,105),xlim=c(-80,50),ylab="",xlab="% change in pact onset hazard",axes=F,pch=19,cex=1.5)
abline(v=0,lty=2)
segments(x0=lower,y0=yrng,x1=upper,y1=yrng,lty=1,lwd=1.5)
labels<-c("Electoral cycle", "Left Gov.", "Unemployment", "Inflation")
labels.sub<-c("1 additional qtr","median seat chg.","median qtrly chg.","median qtrly chg.")
labels.sub.sub<-c("(1/16 CIEP)","(4.3)","(0.2)","(0.3)")
text(x=-50,y=c(100,75,50,25),labels=labels,font=2)
text(x=-50,y=c(95,70,45,20),labels=labels.sub)
text(x=-50,y=c(90,65,40,15),labels=labels.sub.sub,cex=0.80)
marks<-c("-20","-10", "0", "10", "20", "30", "40")
poss2<-c(-20,-10,0,10,20,30,40)
axis(1,at=poss2,labels=marks)

dev.off()



